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ABSTRACT 


A  computationally  feasible  and  strongly  consistent  method  is  considered  for 
estimating  the  coefficients  of 


Vt  +  V,-i+' 


+  y  x 
P  t-p 


u  . 
t 


The  steady-state  solution  is  observed  in  the  presence  of  nonvvhite  noise,  and  the 
system  is  driven  by  a  certain  selected  superposition  of  sinusoids,  u^.  When  the  noise 
is  white,  the  estimates  are  shown  to  possess  an  asymptotic  normal  distribution  with 
covariance  matrix  which  depends  on  the  coefficients  only  via  the  values  of 

P 

,  V  ikcj  ,  2 

lAV  1 

k=0 


evaluated  at  the  input  frequencies. 
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1.  INTRODUCTION 


We  observe  the  steady  state  solution  of  a  p  order  difference  equation  in  the  pre¬ 
sence  of  additive  noise.  The  deterministic  input  sequence  is  at  our  disposal,  up  to  an 
unknown  gain.  Corresponding  to  a  superposition  of  sinusoidal  inputs,  vve  exhibit  an 
estimate  of  the  p  +  1  parameters  having  the  following  properties. 

(i)  On-line  computation  of  the  estimates  can  be  done  recursively  in  the  incoming 
data  by  operations  which  entail  only  first-order  storage  requirements,  (ii)  As  the 
observation  time  increases  to  infinity  the  estimate  is  strongly  consistent  under  very 
weak  conditions  on  the  noise  process,  which  include  stationarity  but  not  the  existence 
of  a  power  spectrum,  (iii)  The  standardized  estimates  tend  to  joint  normality  with  a 
covariance  matrix  which  is  relatively  easy  to  calculate  when  the  noise  is  white.  If, 
in  addition,  the  noise  is  Gaussian,  the  estimate  has  minimal  (generalized)  variance. 

The  initial  step  is  "naive”  least  squares,  i.e.  solving  normal  equations  (and  hence 
(i)).  Simple  instantaneous  nonlinear  operations  are  then  performed  on  the  least  squares 
solution.  The  result  is  mapped  by  a  fixed  matrix  into  the  estimate  of  the  ratio  of 
coefficients  to  gain. 

2.  THE  MODEL  AND  ESTIMATE 


Corresponding  to  a  given  deterministic  sequence  {u^},  let  { x^}  < 


satisfy 


x  =  0.  x  +*-*+0  x  +  Ku  (t  =  . .  .  , -1,  0,  1,  .  . . ) 
t  1  t-1  p  t-p  t 

for  some  input  gain  K  f  0.  We  will  assume  stability,  i.e.  that  the  p  roots  of  the 
characteristic  equation 


are  each  in  modulus  less  than  unity.  We  will  also  assume,  but  only  for  the  sake  of 
simplicity,  that  the  roots  are  distinct. 

We  observe 

yt  =  xt  +  (t=  2 . n,  .  .  . ), 


1 


where  {^}  is  some  zero  mean  stochastic  process,  and  suppose  that  the  regression  se¬ 
quence  was  arbitrarily  initialized  in  the  remote  past.  In  statistical  terminology  we 
are  faced  with  a  particular  nonlinear  regression  problem: 


x  =  F  (K,  0  0  ) 

t  t  1  p 


is  known  up  to  the  values  of  p  + 1  parameters  to  be  estimated. 

The  input.  To  estimate  K  and  9^,  . .  . ,  6  we  distinguish  between  odd  and  even 
values  of  p  +  1,  and  accordingly  take  for  the  input 


1 


+ 


/  cos  at  t 
Lj  J 
j  =  l 


u 

t 


+ 


V  (-1) 

N  cos  a)  t  + 

U  J 


t 


1=1 


*7  2 


if  p  =  2q 


if  p  =  2q+  1. 


The  angular  frequencies 

0  <  a:  1  <  •  •  •  <  at  <  7t 

are  to  be  chosen.  Thus,  the  input  is  l/\/2  when  p  =  0,  giving  the  constant  regression 
K/'s/ 2  for  the  trivial  case  of  no  O' s.  The  input  is  1/^/2  +  (-l)V\/2  when  p=  1, 
l/x/2  +  cos  Ct^t  when  p=  2,  and  so  on. 

We  present  the  estimation  procedure  for  an  even  number  of  parameters.  After 
doing  so,  we  indicate  the  modifications  to  be  made  for  an  odd  number,  as  well  as  for 
the  case  when  the  gain  is  known.  Thus,  in  what  follows, 


p  + 1  =  2q  +  2 


for  some  integer  q  ^  0. 
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Step  1.  Choose  q  distinct  frequencies  interior  to  (0,  tt),  and  define 


h' 


t 


1  .  (-if 

-  ,  costa.t,  sinca.t,  .  . .  , cos co  t,  sina;  t, - 

^2  1  q  \I~2 


for  t=  1,  2,  .  .  . ,  n  with  n  >  p.  Define  matrices 


H 

n 


h  1 
n 


B  =  (H  H'  ) 
n  n  n 


where  prime  denotes  transposition.  Let  be  the  column  vector  of  the  first  n 

observations  ,  and  let 

i  n 


$  =  R  H  y  . 

n  n  n  n 


It  is  known  that  successive  vectors  so  defined  (the  solutions  of  normal  equations  for 
increasing  sample  size)  satisfy  the  difference  equation 


(P 


<p  +  k  (y  -  h 9  <p  ,) 
n-1  n  n  n  n-1 


(n  =  1,  2, . . .  ;  =  0) 


where  the  gain  vectors  can  be  computed  recursively  without  matrix  inversions.  We 
are  not  interested  in  the  details  of  such  computations  -  only  that  they  can  be  done. 

Step  2„  At  each  n,  then,  we  have  a  (p+1)  -  vector  whose  components  (for 
reasons  which  will  become  clear)  we  label  as 


< V 


/ 

n 


1 


From  it  we  compute  a  new  vector 


r 

n 


=  r 

L 


(n)  (n)  » 


'0 


1 


1 


c(n) 

q+1  • 


3 


by  means  of  the  relations 


These  hold  for  j  =0,  1, . . . ,  q,  q+  1  with  bj;n)  =  b(n), 

0  q+1 

Step  3,  Next  form  the  p  +  1  by  p  +  1  (p  =  1,  3 , . . 


=  0. 

. )  matrix 


1  1  1 

1  costa^  cos2ca^ 

0  sinta^  sin2o>^ 


1 

cos  pea  ^ 
sin  pea  ^ 


M  = 


1 

cos  eUq 

cos  2c a 

COS  p CO 

0 

sin  eu^ 

sin  2ea^ 

sinp co 

1 

-1 

+  1 

+1 

This  matrix  is  nonsingular,  and  we  can  compute 


4 


y 

n 


r 


y 


(n) 

0 

(n) 

1 


=  M  ip 


for  n  >  p. 

Step  4.  The  estimate  of  K  based  on  the  first  n  observations  is  1/y 
0^  is  -y^/  y^  (k=  1,  2 . p).  If  the  original  system  is  rewritten 


(n) 
0  ’ 


and  that  of 


Y0xt+Ylxt-1+' 


+  y  x  =  u 
P  t-p  t 


and  it  is  the  y's  which  are  to  be  estimated,  then  the  procedure  terminates  with  Step  3. 
Odd  number  of  parameters.  When 

p  +  1  =  2q  +  1  (q  >  0) 


we  modify  the  preceding  by 

(i)  deleting  the  last  entry  of  h  in  Step  1, 

(ii)  deleting  the  last  entries  of  <p  and  ip  in  Step  2, 
(iii)  deleting  the  last  row  of  M  in  Step  3, 


and  leave  everything  else  unchanged. 

Choice  of  input  frequencies  .The  inversion  of  M  is  obviated  when  we  use  the  par¬ 
ticular  angular  frequencies 


_  27Tj 
P  +  1 


q)- 
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For  p  +  1  =  2q  +  2,  put 


1 

p  +  1 

2 

p+1 


0 


s  = 


0 


2 

p  +  1 


1 

p  +  1 


The  matrix  Sz  M  is  orthogonal, 
Step  3.  We  get  for  k  =  0,  1, . . . ,  p 


so  it  is  any  easy  matter  to  solve  the  equations  in 


/  (n)  (n)  .  . 

i  c  coskca.  +  d.  sinkcu. 

V  i  j  j  J 


(-1) 

+ 


c 


(n) 
q+1  ' 


j=l 

Forp  +  l  =  2q+l  the  last  term  is  dropped. 

Known  input  gain.  When  K  is  known,  we  replace  the  input  u^  by  u  /K  and  every¬ 
where  read  p  for  p  +  1.  In  Step  3  the  first  column  of  M  is  deleted,  and  the  estimate 

of  0  is  the  k^1  component  of  -M  ^  ip  (k  =  1,  2, .  . . ,  p),  i.c.  y^  =  1. 

K  n  u 

3.  THE  ASYMPTOTIC  BEHAVIOR  OF  THE  ESTIMATES 


The  following  two  lemmas  preceding  the  theorems  of  interest  contain  known  re¬ 
sults,  and  are  included  only  to  make  the  report  self-contained.  Proofs  of  the  lemmas 
and  theorems  are  given  in  Sec.  4  and  Sec.  5  respectively. 

Lemma  1.  Let  Qt, . . .  ,  6  be  any  real  numbers  such  that  the  roots 

-  1  p 


of 


1  p 


g(X)  =  XP-01XP_1 


0=0 

P 


6 


are  in  modulus  less  than  unity  and  distinct.  Let  { v^}  be  any  bounded  number  se¬ 
quence.  Then  the  steady  state  solution  of 


is 


x  =  0,  x  .+••-  +  0  x  +  v 
t  1  t-1  p  t-p  t 


x  =  /  B  v 

t  /j  n  t-n 


n=0 


where 


K  “k\" 


k=l 


kP"1 


O'  = 

k 


<W 


Lemma  2.  For  the  particular  input 


v^  =  K  coscat  (K^O,  ca  arbitrary) 


the  conclusion  of  Lemma  1  reads 


where 


and 


x  =  a  cos  cat  +  b  sin  cat 
t 


a  +  b 

b 


.  y  cos  kca 
2  .2  /  ;  k 


k=0 

P 


a  +  b 


-z 


2  .2  /,\Slnk“ 


k=0 


yn  “  k 


\  =  -  T  <k  =  1’2 . p)' 
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Theorem  1.  In  the  notation  of  Sec.  1  let 


O' 


O' 

n 


[K,  0. . 0  } 

1  1  p  1 

r  ,  ^n)  r<n)  i 

lT"*”'  "T"7 ’3^ 


with  p  +  1  even  or  odd  as  the  case  may  be.  Then 


P { lim  0  =(?}  =  ! 
_  n 


for  any  distinct  interior  points  cOj,  .  . .  ,  of  (0, 7r)  and  any  0  mean  stochastic  process 
{£  }  for  which 

o(i/h£) 


as  h  — *  00  for  some  e  >  0. 


Theorem  2.  Let  {|  }  be  a  sequence  of  independent  zero  mean  random  variables 
^2  C  2  +  6 

with  common  variance  cr  ,  and  suppose  sup  £  |£  |  <°°  for  some  6  >  0.  Choose 

t  1 


0  =  o)Q  <  cu1  < 


<  CtK  <  CO 


q+1 


=  IT, 


which  are  otherwise  unrestricted. 

(a)  The  (p  +  l)-vector  \m( y  •y)  has  a  large  sample  normal  distribution  about  the 
origin  with  covariance  matrix,  for  p  +  1  =  2q+  2  even, 

Z  =  2<t2M"1P2m'1' 

P  =diag[p(),p1,p1 . Pq.Pq.Pq+1  ]• 

M  is  given  in  Step  3,  and  the  entries  of  P  are 
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0=0,1 . q,  q+1) 


rte 


ikco . 

=.  J  I2 


k=0 


(b)  The  covariance  matrix  of  the  asymptotic  normal  distribution  of  ^n (0  -  0) 

n 

is 


A'  2  A 


where 


I 


0 

P 


and  I  is  the  px  p  identity.  The  entries  of  P  are  also  expressible  as 

.P,  ikco .  icu 

pi  =  ~r j  I  =  ~T  ls(e  J)  1 2  ( g  =  g  )  » 

1  U  k  K2  p 

k=l 

i.e.  in  terms  of  the  coefficients  0A,  . . . ,  0  or  the  characteristic  roots  A  .  .  .  . ,  A  . 

1  p  1  p 

(c)  For  p  +  1  =  2q  +  l  odd,  we  delete  the  last  row  of  M  and  the  last  entry  of  P 
under  (a). 

Corollary.  For 


co 

j 


( j  =  1,  2,  .  .  . ,  q) 


the  covariance  matrix  2  in  Theorem  2(a)  becomes  a  Toeplitz  matrix  with  entries 


a 


ab 


2a 


2 


(p+1)2 


q 

P0  +  4  ^  Pj2cos(b-a)co.  +  p2+1(-l)b  a 

j=l 


(a,  b  =  0,  1, ....  p). 
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The  generalized  variance  of  the  asymptotic  distribution  of  %/n( y  -y)  is 


da,  r  -  .■ 


2  \',+  ]  2  4  4  2 


^  P  +  1  J  P0/Jl  PqPq+l  ’ 


and  that  of  \n (0  -  0)  is 
n 


/  4,t“  XP+1  i  '->224  242  2 

4t  —  —  <K  P„>  <K  Pj[)  •••  (K  pq)  (K  pq+1)  , 


p+V  K2p  n 


2 

Each  of  the  factors  K  p  depends  only  on  the  0's,  and  is  no  larger  than  4P. 

Confidence  intervals  free  of  unknowns.  We  can  consistently  estimate  p  by 


(n) 

•i  ■ 


«■’)  *  oo 


2  (j  =  0,  1,  . .  . ,  q,  q+1  ;  bj"'  «  b*”^  =  0) 


(n)  (n)  (n)  (n)  (n)  (n)  , 

resulting  from  Step  1.  Letting  be  the  diagonal  matrix  [  p^  ,  p^  ,  p^  . pq  .  Pq  ,  pq+^  ], 

the  vector  of  p+1  =  2q+2  random  variables 


P  M  \  n  (y  -y ) 
n  n 


z  (y-h'<P  ) 

n-l  Lj  t  t  n 

t=l 


has  a  standardized  multivariate  normal  distribution  in  large  samples. 

Again,  it  is  only  necessary  to  delete  all  terms  involving  p  when 

p  +  1  =  2q+  1,  as  well  as  the  last  row  of  M. 

K  known.  The  matrix  2  is  now  px  p  with  2q  +  2  is  redefined  as  p  when  it  is  even. 

Under  Theorem  2(a)  the  first  column  of  M  is  (has  been)  deleted,  and  in  the  definition 

of  p.  we  set  y  =  1  and  y  =  -0  (k  =  1,  2, . .  .  ,p).  Under  Theorem  2(b),  we  set- A 
J  d  k  k 

equal  to  the  pxp  identity  matrix,  and  K  equal  to  unity  in  the  formula  for  p..  Thus,  in 
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the  Corollary,  the  selected  frequencies  are  co  -  2tt j / p.  In  the  formula  for  <7  ,  ,  the 

22  J  ab 

multiplier  becomes  2a  /p  and  a,b  runs  over  1,2,  . .  .  ,  p. 

Remark.  When  the  independent  errors  are  normal,  0  is  the  Maximum  Likeli- 
-  n 

hood  estimate  of  0  for  every  n.  This  is  true  since  the  Least  Squares  estimate  <p  of 
iP  becomes  the  Maximum  Likelihood  estimate,  and  the  Maximum  Likelihood  estimate 
of  any  1-1  vector-valued  function  is  the  function  of  the  Maximum  Likelihood  estimate. 


4.  PROOF  OF  THE  LEMMAS 

t  h 

Proof  of  Lemma  1.  In  the  standard  fashion  we  write  our  p  order  equation  as  a 
first  order  vector  equation: 

x  =  0  x  +  v  e , 
t  t-1  t  1 


where 


Vi 

r  1 

e.  ...  0  • 

1  p-i  , 

0 

p 

n 

1 

4 

n 

Vi  = 

• 

el  = 

• 

e  = 

1 

I 

1 

. 

vp 

n 

1 

1 

0 

Iterating  back  to  t  =  0,  the  solution  for  positive  times  is 


where 


x 

t 


e' 

1 


x 

t 


=  e' 


+ 


t-i 

i_i  n  t-n 


n=0 


J8  =  e'.  0  e 
n  1  1 


We  first  use  the  assumption  of  distinct  roots  to  evaluate  /3^. 
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Throughout  A  denotes  an  indeterminate  complex  scalar.  Define  the  unimodular 
(i.e.  determinant  =  1)  matrix 


E(A)  = 


A 


2 

A 


g'1 


1 


A 


AP'2 


0 


1 


A 

i 


Then 


E(A)  =  E^AJE^X)  ...  E  (X) 


where  E.(A)  is  the  identify  matrix  with  A  added  into  the  (j,j  +  1)  position.  Right- 


multiplication  of  any  p*p  matrix  A  =  [a  ,...,a  ]  by  E.(A)  adds  A  times  the  j 
.  .  .  ,st  „  .  P  J 


th 


column  to  the  j  +  1  .  Consequently, 


AE(A)  -  [  a^,  +  a^A,  a.^  +  a^A  +  a^A2 . ap+an-l*  +  •  •  •  +  AP  ] 

where  the  successive  columns  are  A-polynomials  with  vector  coefficients.  The  in¬ 
verse  of  the  operator, 

E(A)”1  =  E^A)  +  E2(A)  +  -*-  +E  (A), 

exists  independently  of  A.  If  we  replace  A  by  the  characteristic  matrix  of  O  we  find 


( 0  -  A I  )  E(A)  = 


g4(M  k2(a) 


gP-iw 


g  (Al 
^P 


0  J 


=  F(A) 
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wherein 


g.(A)  =  Aj-01Aj'1 . 0  4X  -0.  0  =  1,2 . p) 

and  g^  is  g  in  the  statement  of  the  lemma.  Since  |  E(A)|  =  1,  9  -XI  and  F(A)  are 
rank  equivalent,  and  the  characteristic  polynomial  of  9  is  simply 

|  0-AI  |=  |  F(A)  |  =  (-l)P+1gp(A)  . 


Regarded  as  a  function  of  A,  the  characteristic  matrix  has  full  rank  unless  A  coincides 
with  one  of  the  roots 


of  gp(A)  =  0.  in  this  case,  0  -  A  I  has  rank  p-1.  Thus,  if  A^  is  repeated  there  is 
only  one  nonzero  vector  solution  to  the  homogeneous  equation  (O  -A  l)x  =  0.  Consequent¬ 
ly,  0  is  equivalent  (more  properly,  collinear)  to  the  diagonal  matrix  A  of  its  eigen¬ 
values  if  and  only  if  the  eigenvalues  are  distinct.  This  is  our  assumption. 

Since  0  is  not  normal,  i.  e.  does  not  commute  with  its  transpose,  two  sets  of 
eigenvectors  are  involved.  To  each  A  belonging  to  the  spectrum  A^,  . . . ,  A^  there 

are  distinct  nonzero  vectors  l  -  l  and  r  =  r  ,  unique  up  to  a  multiplicative  constant, 

A  A 

satisfying 


9r  =  Ar 


We  use  the  "left-right"  nomenclature  relative  to  0.  Let 


L  =  [  l  , . . . ,  l  ]  R  =  [  r . r  ] 

1  1  D  1  P 


be  the  matrices  of  these  columns  generated  as  A  assumes  the  p  spectral  values. 
Then 

0R  =  0  [  r  ,  . .  .  ,  r  ]  =  [Or  . Or  ]  =  [X  r  , . . . ,  A  r  ] 

1  1  pJ  1  1  p  L  1  1  P  P 

=  RA 
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and 


O'L  =  L  A  or  L,Q  =  AL' 


Thus,  (without  restriction  on  the  eigenvalues) 


A(L'R)  =  L'QR  =  (L'R)  A 

where  L'R  has  the  inner  product  -tfr.  for  ij  entry.  In  terms  of  elements,  the 
statement  that  L'R  commutes  with  A  is 


(A.  -  A.)  V.  r .  =  0  . 
i  J  i  J 

Thus,  for  distinct  eigenvalues,  L.  is  perpendicular  to  r^  for  all  i  f  j  (and  the  converse 
is  also  true).  The  two  vector  sets  are  called  bi-orthogonal,  and 


L'R  =  D  = 


11 


22 


PP 


d  .  =  L'r. 
JJ  j  J 


We  then  have,  since  I  D  |  J-  0, 


G  =  RAR 


-1 


_  i 


=  (RD  2  )  A(RD  2)'  , 


and  the  proper  normalization  for  t.  and  r.  is  l/'v/d” .  . 

J  J  JJ 

We  now  compute  L,  R  and  D  as  functions  of  the  eigenvalues.  We  have 
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(O' -MM  =  E(A.)'_1  F(A.)'  t.  . 

J  J  J  J  J 

By  definition  -f .  makes  the  left  side  the  0  vector,  and  hence  F(A.)'l.  = 
J  J  J 

equations  are 


-g^j)  1 

J  1 

-g2(\j)  1  1 

1 

gi(V 

g2<V 

1 

>_  i(V '  1 

• 

0  1  0  ...  “0 

because  g  (A)  =  0  for  A  =  A.. 


The  vector  equation  Or.  =  A.r.  is 

J  J  J 


e  0  ...  0  1  B 

1  2  p-1  1  P 

AP_1 

j 

1 

1 _ 

1  ~  -  -  - 

1  1 

AP'2 

J 

= 

AP_1 

J 

.  1 

1 

1  1  0 

1 

1  _ 

1 

_ 1 

The  first  of  these  says  A.  is  a  root,  and  the  balance  are  redundancies. 


i  ~ 

J 

xp"4 

J 

l.  = 

gi(V 

r  = 

AP‘2 

j 

J 

j 

1 

0.  These 


Thus, 
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R  is  a  Vandermonde  matrix  with 


R 


M-ir1  it  w 

i  <  j 


To  compute  the  normalization 


P-1 


d..-t!r  =XP‘‘  +  >  gl(X.»rfc-1 
JJ  J  J  J  L  k  j  j 


k=  1 


we  derive  a  A-identity.  If  we  multiply 


h w  =  /-zvl'‘ 

i=l 

by  A^  ^  and  sum  over  k  we  have 


p-1 


p-1  k 


g|twxp-|c-1  =  (p-i rf-'S  V» 


k=l 


6 .  A 

i—J  L-!  i 
k=l  i=l 


=  pA 


p-i 

p-i_xp-i_  y 

x_ 

k=l 


Ok(p-k)A 


p-k- 1 


Moving  A^  1  to  the  left  side  we  obtain 

»p-1  +  lWk-1  ■  1 

k=l 


p-1 

Ap.y„  >p-k+cl 
L  k 

k=l 


We  can  set  the  arbitrary  constant  equal  to  -0  .  There  results 
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xP-1+V^<»xP-k-*  ^gp(X) 


k=l 

p  p 


■^n»-v  - 1  n  <*-»,> 

k=l  i  f  k 


i=l 


In  particular,  for  A  =  A.  , 


d„  =  TT  (  X  -  A  )  . 

JJ  .j .  J  1 
1  *  J 


To  find  the  dependence  of  j3  on  the  eigenvalues  we  use  L'O  =  AL ' , 


A>' 
1  1 


L'Qn  =  An  L'  = 


*2*2 


A*V 

P  P 


Since  the  leading  element  of  every  t  is  1, 


L'0  e  = 
1 


-1  -1 

But  L'R  =  D  implies  L'  =  RD  ,  so 


We  have 
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0  el  = 


RD 


-1 


An 

LP  J 


Since  e'R  is  the  first  row  of  R, 
1 


=  K' 1 


A^"  1  ]  D_1 
P 


and  therefore 


p  ap-i 


= 


J 


d..  J 
i-1  JJ 


A"  = 


a  A 
J  j 


i=l 


The  vector  ot  satisfies 


L  ot  =  e 


In 


particular,  ot  ^  +•  •  •  +  ot  =  1,  as  should  be  since  /3^  must  be  unity. 
We  now  use  the  assumption  that 


max  i  A.  |  <  1 


1  <i  < 


J  -P 


J 


where 

norm 


denotes  modulus.  For  any  p*q  matrix  A  with  elements  a_ 


,  define  the 
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which  generalizes  length.  Then 

e'l  0,1  ><n  |  s  ||  e’  j|  ||  enx0ll 

=  i|0n||  l|x0l|  . 

But 

II  enii  ||  R  ||  ||  An||  ||  R- *||  =  const.  J  X^n  +  ---+XV 

1  P 


— ■  0  as  n  —  °°  . 


Translating  the  origin  back  to  t  =  -  °°,  this  establishes  the  formula  asserted  in  Lemma  1. 

Proof  of  Lemma  2.  For  v  =  Kcoswt  the  steady  state  solution  is,  according  to 
Lemma  1, 

=  (Ka')coswt  +  (Kb')sinwt 

where 


a'  = 


B  cos  nw 

n 


=  a/K 


n==0 


b1  =  /  B  sinnw  =  a/K  . 
Lj  n 

n^O 


Setting 


A  =  1  -  N  0kcos  kw 

k=l 


B  =  >  fl^sinkw  , 

k=l 
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the  assertion  to  be  proved  is  clearly  equivalent  to 


2  2 
a '  +  b' 


=  A 


b' 


2  2 
a'  +  b' 


=  B  . 


This  in  turn  holds  if 


(A  +  iB)(a'-ib')  =  1, 


since  then 


A  +  iB  = 


1 


+  ib' 


a ' -  ib' 


+  i- 


+  b' 


+  b' 


Putting 


icu 

z  =  e 


we  have 

P 

A  +  iB  =  1  -  /  8  z  ^  =  z  (z) 

k  V 


and 


a'  -ib' 


n>0 


P  P 


n^O  k=l  k=l  nS:0 


But  |  A  /z|  s 
k 


<  1  by  hypothesis  for  all  k, 


a '  -  ib ' 


P 

=Y 

£ 


a 

k 


k=l 


with  a  = 
k 


20 


It  thus  suffices  to  prove 


'  Q' 

/  \  \  k  _  p-1 

g  (z)  - —  =  z 

p  Lj  z-X, 

k=l  k 


Since  g  (z)  =  (z-X.),  the  left  side  is 


i  =  l 


L(z) 


k=l  j^k 


<*-v  -ir  n 


z-  X. 

_ I_ 

X,  -  X. 

k=l  J  ^k  k  j 

Ak_1  Lk(7)  * 
k=l 


Each  of  the  summands  is  a  (Lagrange)  polynomial  of  degree  p-1  in  z,  so  we  must 
have 


L(z)  =  c^t} 


p-2 

+  c2zP  +  ' 


+  c 


The  coefficients  are  determined  by  the  values  of  L  at  any  p  distinct  points.  We  of 
course  use  the  characteristic  roots,  since 


Thus 


L.  (X.)  = 
k  1 


1  if  k  =  i 


0  if  Mi 


L(X.) 


This  shows  c  =  1  and 
1 


c„  = 


XP  (i=  1,  2, . . . ,  p)  . 

c  =  0,  i.  e.  L(z)  =  z^  \ 

P 
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5.  PROOF  OF  THE  THEOREMS 


Proof  of  Theorem  1.  It  clearly  suffices  to  consider  either  an  odd  or  an  even 
number  of  parameters,  and  we  choose  the  latter.  Thus,  throughout  the  proof, 


p  +  1  =  2q+  2  (q  2  0). 


Furthermore,  we  set 


"o  =  0 


co  =7 r . 

q+1 


According  to  Lemma  2,  then,  the  output  which  results  from  the  input 

q 

i  v  1 

u  =  -  COS  t  +  ;  COSW.t  + - COSCO  ,t 

c  -Jl  0  L  J  ^  q+1 

J=1 


IS 


x  = — -(ancosw  t  +b.sinw  t)  +  /  (a. cos w.t  +  b.sin w.t) 

t  i-  o  o  o  n  i_j  j  j  j  j 


*7  2 


j=l 

1 

+ - (a  ,  .  cosca  ,  .t  +  b  sinca  .  .t  ) 

q+1  q+1  q+1  q+1 


=  h' >  , 


where  h  is  as  in  Step  1  and 


cp*  [a_,  a  ,  b.,.,,,a  ,b  ,a  ] 
L  0  1  1  q  q  q+1 


Furthermore, 


c  = 

j 


J 


2  J- 

a.  +  b. 
J  J 


=  )  v,  cos  kca. 

U  k  j 

k=0 


j  =  0,  1,  .  .  .  ,  q,  q+1 


d.  s 
J 


2  2 
a.  +  b. 

J  J 


7,  sinkta. 
k  J 


k=0 
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which  imply  h^  =  b  ^  =  0.  Putting 


*'  =  [co'crlli . W  V*1 

y'  =  [r0»r1»---.r  ] » 


where  yn  =  1/K  and  y  =  -  0  / K  (k=l,  2, ....  p),  these  equations  are 

v  t  K  K 


ip  =  My 


with  M  as  given  in  Step  3. 

Let  £  be  the  columr 
n 

model  which  is  linear  in  the  components  of  cp, 


Let  £  be  the  column  vector  of  the  first  n  noise  realizations.  We  have  i 
n 


y  =  H'  <p  +  | 
n  n  n 


The  estimate  resulting  from  Step  1  can  be  written  <p  =  <p  +  B  H  £  ,  or 

n  n  n  n 


A  (<p  -  <p)  =  H  £ 
n  n  n  n 


(A  =B  ). 
n  n 


For  the  p+1  by  p+1  matrix  A  =  >  h  h'  we  have 

n  Lj  t  t 

t=l 


data 
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2  2  n 
-A  =-  2 
n  n  nt=1 


i  1  1  . 

•>  - cosco  t  - sinco t 

■J2  1  -12  1 


n/  2 


■cos  co^t 


cos  to  t  costo  tsinto  t  costo  tcosto  r 
1  11  12 


sin  co  t  sinto  tcosto  t 

1  12 


cos  co  t 


(symmetric) 


sinco  t 

41  q  *n/2 


.  cosco  tsinco  t  cosco  t 
1  q  1 


(-1) 


.  sinco  tsinco  t  sinco 

1  q 


(-D 

1  41 


.  2  (-1) 
sin  co  t  sinco  l - 

q  q  42 


From  the  identities 


n 

'S  cos  2At 

L-J 

t  =  l 


sin  nA 
sin  A 


cos(n  +  1)A 


n 

V  .  sinnA  .  y 

)  sin2At  =  — : - sin(n  +  l)A 

l_j  sin  A 

t=l 


(A  f  a  multiple  of  i r) 


it  follows  for  n  -*■  00  that 
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2  \ 

—  cos  X  .t  cos  A,,t 
n  .  j  1  2 


t=l 


)  =  0(  — )  for  all  0  ^  7r  with  A .  ^  A„ 

n  12  12 


—  sin  A  t  sin  At 
n  l_,  1  2 

t=l 


n 

2  V  2  t 

—  ;  cos  A 
n  Lj 

t=l 


n 

2  V  .  2  t 

—  /  sin  A 
n 


}  =  1  +  0( —  )  for  all  0  <  A  <  it 

n 


t=l 


n 

2  V 

—  '  cosA,tsinA0t 
n  i_j  1  2 

t=l 


°(i) 


for  all  0  <\\  ^  7r 

1  2 


Since  . . .  ,  cj  were  chosen  as  interior  points  of  (0,  7r),  we  therefore  have 


I  +  -E 
n  n 


for  some  matrix  E  whose  elements  remain  uniformly  bounded  as  n  —  00 .  (To  get  the 
n  _ 

identity  matrix  in  the  limit  was  our  reason  for  using  the  amplitude  i/\f2  for  the  first 
and  last  input.)  Thus, 


(I  +-E  )  (0  -cp)=- H  r 

n  n  n  n  n  n 


and  the  estimate  cp  will  converge  to  0  in  the  same  probabilistic  sense  that  the  vector 
n 
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—  H  £  converges  to  the  zero  vector. 

n  n  n 

2  — 

The  p+1  components  of  —  H  £  are  arithmetic  means 

n  n  n 


4i  \ 


n  t 

t=l 


n 

2  V 

—  cosca  tf 
n  Zj  I  t 

t=l 

n 


_  ,  sinca  t£ 

n  i  t 

t=l 


COSO)  tf 

n  q  t 

t=l 

n 


sin  to  t  £ 

n  l_i  q  t 

t=l 


t=l 


each  of  which  has  zero  expectation.  For  any  average  of  the  form 

n 


S  =  —  7  cos  cat  £ 

n  n  /_j  t 

t=l 
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we  have 


n 

c  1  V 

(  S  cosconf  =-  cos  cot  cos  con  a 

n  n  n  ,  ,  n-t 


t”  1 


By  hypothesis,  therefore. 

]  rCS  cos  co 


n- 1 

iVi 


1 


a  =  Of - ) 

n  n  '  t 1  € 

t=o  n 


It  follows  from  a  LLN  for  dependent  random  variables  that  0  with  probability 

one  as  n  —  00  (anti,  incidentally,  also  in  mean  square). 

Since  strong  convergence  is  preserved  through  continuous  transformations,  we 

successively  conclude  from  to  —  that  d>  —  ib  ,  y  —  y,  and  the  conclusion  of  the 

n  n  n 

lemma. 

Proof  of  Theorem  2.  We  continue  to  take  p+1  even,  and  note  the  modifications 
of  the  results  needed  for  p+1  odd  when  we  are  done.  The  asymptotic  distribution  at 

__  2  _ k 

\in(<p  -  cp)  is  the  same  as  that  of  — —  H  £  ,  which  we  express  in  symbols  as 

n  /—  n  n 

V  n 

\  n(c/?  -  < v )  ~  - H  t; 

n  r  n  n 

Vn 


Under  the  assumptions  of  the  theorem  the  latter,  by  Liapounov's  CLT.  tends  to  be 
distributed  in  large  samples  as  ap+1  dimensional  normal  random  variable  with  zero 
mean  vector  and  covariance  matrix  equal  to 


lim  (' 
n 


( 


2  —  \ 
-H  | 


n/  n 


—  n  n 


H  l 


n  n 


t 


According  to  the  asymptotic  orthogonality  relations  used  in  the  proof  of  Theorem  1, 
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n  n 


C-4  V  V  .  t  4<r  V 

(  —  /  /  COS  Ca  t  £  COS  Cl)  s|  = -  >  COS  Cl)  t  COS  Cl)  t 

n  L  L  i  t  j  s  n  l_j  i  j 


2  a2  <5. 


ij 


t=l  s=l 


t=l 


as  n  —  00  for  all  i,  j  =  0,  1,  .  .  . ,  q,  q+1,  provided  i  and  j  are  not  both  0  nor  both  q+1. 

The  same  limit  holds  when  both  cosines  are  replaced  by  sines.  For  i  =  j=0  (i.  e.  ca^=0) 

and  i=j  =  q+l  (i.e.  ca  =  it)  we  have 

q+1 


2  V  ft  {  .  £lf  f  *2 


t=l  s=l 


t=l  S=1 


for  every  n.  Finally, 
n  n 


e- 

n 


V 


cos  ca  t  |  sinca.s  £ 

U  it  J  s 

t=l  s=l 


2  n 
4  a  V 

-  /  cosca.t  sinca.t 

n  L  i  J 

t=l 


without  restriction  on  i,  j  =0,  1 . q,  q+1.  Thus, 

\  n (<p  -  <p)  ~  N(0,  2a2 1 ) 

n 

2 

where  I  is  the  p  +  1  by  p+1  identity.  This  checks  with  the  known  result  for  N(0,a  ) 

2 

errors  ;viz.,  that  <p  is  N(<p,  a  B  )  for  every  n.  (Recall  B  ~  2/n. ) 

n  n  n 

Since  the  components  of  ip  are  functions  of  those  of  <p^  we  have,  by  the  ’’delta 
method,  ” 

^n(ipn~ip)  ~  N(0,  2<t  P  ) 

2-12  - 1 ' 

N/n(y  -y)  ~  N(0, 2a  M  P  M  ) 

n 

where 

2  a^ 

P  =  D'D  and  D  =  — 

ocp 
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is  the  p+  1  by  p+1  matrix  of  partial  derivatives  computed  at  the  true  parameter 
values  (induced  by  K  and  the  0's).  We  have 


0c. 

0a. 


Od 


■*-  =  — =  6  cv 

0b.  ij  j 

l  i  J 


■>  •> 
b  -a 

O'  =  J  -  -1  ■ 

j  «*j2+b?2 


0d.  0c. 


— ^  = - *-  ~63 

0a.  0b.  ij  j 

i  l 


2a.  b 

8  =  - U - 

Pi  2  2  2 

(a  +  b  ) 


for  ail  i,  j  =  0,  1 . q.  q+1,  provided  we  continue  to  take  b(J  =  b^+  =  0  and  ignore 

the  variables  not  involved.  We  see  that 


D  = 


O' 


01 


q+1 


•  a 

P.l 

i 

l 

A.  = 

1 

cv_ 

i  i 


Since  the  inner  product  between  any  two  different  columns  of  D  vanishes,  we  have 


2  2  2  2 
P  =  diag[p0>  pj ,  P1 


VVVd  wo= Vi=  0> 


The  typical  entry  is 
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1 


2  2  2  ,2  2 

2  _  (b  -a  )  4a  b 

P  2  2  4  +  2  24 


(a  +b  ) 


2  2  2' 
(a“  +  b  )'  (a  +  b  ) 


But  from  Lemma  2 


2  u2 

a  +b 


2  ,2  ,227 

a  +  b  (a  +b  ) 


P  2  P  2 

y^sinkcu 


i  y,  cos  kca  + 

\  Lu  k  J 


k=0 


/ 

k=0 


In  other  words, 


Pj=lZ,V 

k=0 


ikca.  o 

J  I 


for  j  =0,  1,  . . .  ,  q,  q+4  and  0  =  ca  <  w  <  •  •  •  <  ca  <  ca  =  ir. 
J  1  0  1  q  q+1 


Finally,  we  apply  the  "delta  method"  once  again  to  get  the  large  sample  distribu- 

1  yl  yn 

:s  of  K  =  —  ,  6= - , . . . ,  0  - - .  We  have 


tion  of  the  estimates 


To  1  yo 


y, 


where 


si n(0  -0)  ~  N(0,  2(j2A' M_1  A  ) 

n 


A  = 


9(k,  $4,...,e  ) 

1  p 

8(T0'T1 . V 


=  -K 


where  I  is  the  p*p  identity. 


K 

0 


01  P2 


,  e 
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We  can  give  an  alternate  expression  for  the  elements  of  P  in  terms  of  the 
characteristic  polynomial 


g(z)  =  zP 


1 


Dividing  through  by  K  we  have 

k=0 


ICC 


Putting  z  =  e  and  taking  squared  modulii  gives 


K 


,5lg(e1“')|2  =  |\  Yk. 

k=0 


ikca ,  2 


We  therefore  have 


1  |  .  \  .  ikca ,  2  1  |  ica  ■  2 

pj=  — 1‘- Zske  I  =  -Ir(c  >i 

J  K  ,  K 

k=l 


•  ico  |  tt  j  iuo  j 

One  can  hound  these  quantities  using  the  fact  that  |g(e  )|  =  [I  |e  -  X ^ |  satisfies 

k=l 


nn-iM> s  ig<e‘")i s  n  <‘+hki) 


k=l 


k=l 


For  an  odd  number  of  parameters, 

p  +  1  =  2q  +  1  (q  >  0) 

we  simply  delete  the  last  row  of  M  and  the  last  entry  of  P. 

Proof  of  the  Corollary.  From  the  identities  used  in  the  proof  of  Theorem  1  for 
the  sum  of  the  first  n  cosines  and  sines  we  have,  with 
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<j  =  0’1 . q’ t|+l  =  )  ■ 


relations 


P  P+1 

\  cosu)#k  cos  u).k  =  cos  ca  k  cos  ca  k  +  1  -  co  s  a;,  (p+1)  cos ca (p+1) 
is  1  J  Is  i  j  i  j 

k“0  k=l 


=  Etl 


0  6.. +  1-1 

2  ij 


for  all  i,  j,  provided  both  are  not  0  or  q+1.  From  similar  ones  involving  sines,  and 
cosines  and  sines,  we  see  that 

MM'  =  S"1 

where  S  was  defined  in  Sec.  2.  It  is  now  a  simple  matter  to  solve  My  =  ip  .  Setting 

2 

M  =  S2  M  , 

we  have  M^ M^  =  I,  so  that 

i 

y  =  M^ S2  ip  =  (SM)'  ip  . 


Premultiplication  by  the  diagonal  matrix  S  multiplies  the  rows  of  M.  The  scalar  y 

K 

is  the  inner  product  between  row  k  of  (SM)',  i.e.  column  k  of  SM,  and  ip  (k  =0,  1,  . .  . ,  p). 
Thus 


\  = 


+  • 


p+1  p+1 


(-1) 

(c  cos  ko;  +  d  sinko;)+- — -  c 

j  j  j  j  P+i  q+i 


j=i 


For  p  +  1  =  2q  +  l,  the  last  term  is  deleted. 

The  covariance  matrix  2  under  Theorem  2(a)  becomes 
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1  2  !  2 
— -  -  =  MJS“P  S“MQ  =  M'SP  SM  =  (PSM) '  (PSM) 


The  entry  labeled  a,  b  (a,  b  =  0,  1,  . . .  ,  p)  is  the  inner  product  between  the  corres¬ 
ponding  columns  of  PSM.  Thus 


A  =  fJSL~)2  +  ( 

ab  ^  p  + 1  \  P  + 


2cr 


2P*  \  2 


—  \  cos  ac«;  cos  bca  4*  sin  a cu  sin  bo;  1  +  *  •  • 

1  1  1  1  1 J 


2P  .  \  2 


•••  +  — r  cos  acc  cosbca  +  sinacu  sinbcc  1  +  (-l),+*'> 

V  P+1  >  Q  q  Q  q1  \p+ly' 


which  is  the  asserted  formula.  Again,  if  p  +  i  =  2q  +  l,  the  last  term  is  deleted. 

2  J  2  ! 

The  eigenvalues  of  2  are  clearly  those  of  2<r  S2P  S2  ,  viz. 


r  2  2  2  2  2  2  , 

p  +  1  ‘P0’  2pl  ’  2pl . 2pq’  2pq  ’  Pq+1  ^  * 


The  determinant  is  their  product: 


det 


2a2  \p+'  2 


P  + 


2  2 


i)  »0-  2p4  •  2Pl---  2pq-  2pq-pq+i 


2  v  p+1  „ 

4ct  \ H  2  4 


p  +  1 


4  2 


P0P1  •••  Pq  Pq+1 


Since 


deti  =  (-K)P+’KM-1)P+,KP+2 


we  have 


z  .2  N  p+1  . 

det  (A'  SA)  =  j  (  p 


2  2  2  4 

2p  '  P0^  Pl^ 


2  4  2  2 

<k  <k  pq+,> 


2  ♦ 

The  quantities  K  p.  depend  only  on  the  0's.  Using  the  formula  in  terms  of  the  A's  we 


2  ,p 
see  that  K  p.  <  4  . 
J 
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